Unraveling quantum dissipation in the frequency domain 
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We present a quantum Monte Carlo method for solving the evolution of an open quantum system. 
In our approach, the density operator evolution is unraveled in the frequency domain. Significant 
advantages of this approach arise when the frequency of each dissipative event conveys information 
about the state of the system. 

PACS: 42.50.Lc 



Irreversibility may be incorporated in quantum theory 
by coupling a system to a Markoffian reservoir and trac- 
ing over the reservoir (corresponding to averaging over 
unobserved quantities) to give a description of the sys- 
tem evolution by a reduced density operator equation. 
Recently, a number of theoretical methods have been de- 
veloped which do not perform this trace, but consider 
instead a single trial (quantum trajectory) in which the 
reservoir is continuously monitored Jl|-|^] . Each such trial 
is conditional on a sequence of times t\, t2, • ■ ■, for the 
dissipation events where each t n may be in general associ- 
ated with a decay channel 7„. For example, in the case of 
spontaneous emission from an atom, 7„ would identify a 
unique polarization and direction for the photon. Apart 
from providing valuable insight into the underlying quan- 
tum dynamics, there are significant numerical advantages 
for evolving wave functions rather than reduced density 
operators, and consequently these methods have already 
received widespread application. 

The decomposition of the density operator evolution to 
form a parallel set of quantum trajectories is not unique 
since there are always degrees of freedom associated with 
the quantum measurement basis used to record the ex- 
citations of the reservoir. Significantly, both the insight 
one is able to gain into the dynamics of the system as 
well as the efficiency of the resulting numerical algorithm 
can be very sensitive to this choice. In this letter, we de- 
rive the theory of quantum trajectories (closely related to 
quantum Monte Carlo simulations) in which the unravel- 
ing is done in the frequency domain rather than perform- 
ing the decomposition in time. Consequently, the charac- 
teristic features of previous approaches such as quantum 
jumps of the state and quantum state diffusion are not 
present. Instead we find the observables of the system 
evolve according to a continuous evolution. 

The essential idea is to replace the observed decay 
times t n by frequencies u n . Each quantum trajectory 
is then conditional on a particular record of the reservoir 
state u>i, u>2, . . ., produced by a fictitious measurement 
device of the type illustrated in Fig. H. The output of 
the quantum system is sorted by a cascaded array of fil- 
ters which allow only a particular frequency component 
of the field to pass onto each detector. The unavoidable 



consequence of the filters is that associating a frequency 
with each dissipative event requires that knowledge of 
the precise time at which each decay occurred is lost. If 
we consider a time interval t £ [0, r], in order to form 
a complete reservoir description according to the Fourier 
sampling theorem, each u>i must be chosen from a dis- 
crete but infinite set of frequencies, spaced 2tt/t apart. 
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FIG. 1. The frequency measurement detector. Each output 
channel of the system is sorted by a cascaded array of filters. 

As we will now show, these novel quantum trajectories 
are straight forward to derive. In the case of the measure- 
ment record corresponding to no decays, the choice of the 
time domain or the frequency domain is irrelevant. We 
let |0_r) be the vacuum state in the space of the reservoir 
and denote the state of both and system and reservoir 
by | The quantum trajectory corresponding to no 

decays = (Or \ ®{t)) evolves according to 
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dt ih 

in which the non-Hermitian Hamiltonian H s is related 
to the Hamiltonian for the isolated system -ff sys by 

ih \ -> 
Y 



H e g — H sys — 



(2) 



The operators a 7 (called jump operators in the quantum 
Monte Carlo approach) act in the system space and in- 
duce the corresponding change in the system state when 
a decay into channel 7 occurs. 
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To treat one decay, we introduce the filter operator for 
a single frequency u>i in channel 71 || 

r wi (i) = 4= f e-^^dF^s) (3) 
V r Jo 

where the field operator df 7l (s) acts in the reservoir 
space and annihilates one excitation in the mode 71 in 
the time interval s to s + ds. The quantum trajectory is 
found by applying r Wl (t) to and then projecting 

the filtered state onto the vacuum to enforce one decay 

|^i(*)> = (0fl|r Wl (i)|*(t)>. (4) 

The equation of motion for this is coupled to the evolu- 
tion of \ip(t)) by (see Eq. (153) of Ref. §) 

^|Vw(*)> = ^=\m) + ^(ffeff + M^M*)) (5) 

For two or more decays, one may proceed in different 
directions depending on whether or not time ordering is 
imposed on the dissipative events. We first consider the 
rules for deriving the trajectory in the case of unordered 
measurements (u) . This is the situation applicable for the 
device shown in Fig. |l| where the order of the frequencies 
in the record list plays no role. The associated quantum 
trajectory is defined as 

|v&l.,„„(*)> = (OiiM*) • ••»■«„ (*)|*(*)) ( 6 ) 

which evolves according to 




+ -(H cS + hY / u p )\^l.^(t)) (7) 
P =i 

The first summation couples the trajectory to all n tra- 
jectories which exclude one of the frequencies in the list 
by the associated jump operator for that decay. Note that 
Eq. (||) is a special case of Eq. (|t]) with n = 1. Although 
the evolution of any trajectory can be found by iterating 
Eq. (0) back to Eq. (|j) , the number of coupled equations 
which must be solved grows as 2". Consequently it is 
difficult to treat long time intervals r in which a large 
number of decays may occur in this way. 

The scaling is more favorable in the case of time or- 
dered decays (o) which we now consider. In this case we 
impose the constraint that the w\ decay occurs before the 
L02 decay, which occurs before the W3 decay, and so on. 
Since this corresponds to a different physical measure- 
ment (e.g. an atomic cascade), the quantum trajectories 
are distinct from those of the unordered case. The cor- 
responding nested filter operators are defined recursively 
starting from Eq. (||) by 




The trajectories in this case are 

|^ ) ,..., Wn (*)) = (0fl|r Wl ,..., Wn (*)|*(t)) (8) 
and evolve according to 

||^° 1 ) ,..^(*)>=^I^ 1 ) ,.. J .„_ 1 (*)) 

1 n 
+ - K (H cS + hY / " P )H±..,uJt)) (9) 

P =i 

The difference between this and Eq. (]?]) is that here the 
evolution of each trajectory is coupled only to the trajec- 
tory which excludes the last measurement in the list via 
the jump operator for that decay. The number of cou- 
pled equations which must be solved is equal to n+1. In 
Fig. H we illustrate the solution of this for a few sample 
trajectories for the case of resonance fluorescence from 
a two-level atom. For this case H sys — TiVl{<j + + cr~)/2 
where a + and a~ are the usual raising and lowering op- 
erators and fl is the Rabi frequency. The single jump 
operator ^/Ta~ defines the spontaneous emission rate T. 
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FIG. 2. Sample quantum trajectories |'09r(£)) (dashed), 

|V'o.3r(t)) (dash-dot), and | V'3r,o.5r,-o.3r (t)) (solid), for res- 
onance fluorescence from a two-level atom with 57 = 6F. 

This recipe allows us to calculate arbitrary trajectories 
in both cases, but it remains to show that it represents a 
correct physical unraveling of the reduced density oper- 
ator equation. This is most easily done with the identity 

£|VW..,"n(*)) = ^JVW.-^-lW) (10) 

which may be proved for both unordered and ordered 
trajectories using Eq. (^|) or Eq. (||) to expand the left 
hand side, evaluating the u> n summation to give the Dirac 
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delta function, and applying Eq. (81) of Ref. ||. The 
reduced density operator is constructed by tracing over 
all possible measurement records 



oo 



n— uj\ 



(ii) 



where pu x ,...,u n is 

( (\^ui,...,w n (t))(i>^i,...,u n (t)\)/n\ unordered 



ordered 



(12) 



The n! arises from the trivial permutations of n frequen- 
cies which correspond to the same trajectory. Differenti- 
ating Eq. Jll| ) with respect to time, and using Eq. JTc| ) 
along with either Eq. (0) or Eq. (||) gives in both cases 
the same equation for p 

^=^[#s ys ,p] + ^ ^(2<z 7 pa 7 - a 7 a 7 p - pa^) (13) 
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which is the quantum master equation |6) . This is the key 
result, and it implies that the frequency domain method 
we have presented is a correct unraveling of the reduced 
density operator at all times. 

As is the case in the time domain, the ensemble formed 
by a large but finite set of trajectories, will approximate 
the result of the complete trace and yet involve the evolu- 
tion of wave functions rather than density operators. We 
now outline the simulation procedure we have adopted 
focusing on the case of ordered trajectories. Specifically 
we would like to calculate the evolution of an arbitrary 
observable of the system (O) = Tr{Op} using frequency 
unraveling. The procedure is as follows 

1. Select a system state IV'(O)) from the initial statis- 
tical mixture in p (trivial if p is a pure state) and 
calculate the zero trajectory |^>(t)) from Eq. (Q). 

2. Construct the one decay trajectories cou- 
pled to this with u>i taking on all values from a 
discrete set, {2-kp/t} for p integer. 

3. Using a random number, select a value for oj\ 
weighted by the normalization of the trajectories 
at time r, i.e. from the probability distribution 

PK)= (^.-^W k ll .^(r)) 

£ a , n (Vw..,"»( T ) I Vw..,w„(r)) 

with n = 1 here for the first decay. 

4. Construct \ipoji,u 2 {t)) with u>i fixed and with u>2 
varying over all values given in the discrete set. Se- 
lect a value for ll>2 using Eq. ( |l4| ) with n = 2. 

5. Continue, selecting a frequency for each decay, until 
a predetermined maximum number is reached. 



6. Since we perform a partial evaluation of the series 
in Eq. (pT|), the estimate for (O) is a weighted sum 
over those trajectories which are calculated 

E P t^<^ ) ,..,« B (*)|0|^,.„, tl>B (t)> (15) 

traj 

where Ptraj is the probability of the trajectory being 
evaluated in the algorithm. This probability P tra j is 
unity for the zero and all the one decay trajectories 
which are always calculated, equal to P(wi) for the 
two decay trajectories, P(wi)P(w2) for the three 
decay trajectories, and so on. 

7. Start again from the beginning and average the re- 
sults from many such trials to form the ensemble. 

A simple check of the numerical implementation is to 
establish that Tr{p} = 1 (for O = 1) at all time. In 
Fig. |3| we apply this approach to form the ensemble time 
evolution for the case of resonance fluorescence. 
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FIG. 3. Monte Carlo simulation with 2500 trials for reso- 
nance fluorescence from a two-level atom with Q = 6F. 

Since our method is based on frequency measurement, 
the calculation of spectra arises naturally by bining the 
frequencies of the trajectories in the simulation. For ex- 
ample, we outline now the procedure for using time or- 
dered trajectories to calculate the fluctuation spectrum 

S 7 (w)=lim f f e-"( s - s ' ) (*|dP t (s)dP 7 (s / )|*) 
= lim(*|ri(t)ro,(t)|*> (16) 

which is the rate of decay of frequency ui in channel 7. 
We define a partially ordered trajectory (p) 

l^L^-Jt)) = (OflM*)^,..,^)!^)) (17) 

which evolves according to 



3 



^|^t.,. n _ I W> + ^|^ 1 ) ,.. > ^- 1 ;.(*)) 

+ ^(^cff + hJ2^P + ^) |^ ) ,...^n i <-(*)) (18) 
p=l 

The calculation proceeds in exactly the same way as pre- 
viously outlined except for step 6 which becomes 

6. When calculating | ...,o>„ (t)) with uj n varying 
over all values from the discrete set, evaluate also 
\ipui,...,ui n -t;u n (t)}- The estimate for 5 7 (w n ) is 

E P t^W,...,« B _ I ^ B (T) I ^t.., Wn _ i; .„(T)) 
traj 

where Ptraj is identical to the previous case. 

This computes a transient spectrum unless the steady 
state density operator is used for the initial condition. 
Since r is finite, it should replace the upper limit of both 
integrals in Eq. ( |l6| ) which is equivalent to simulating 
the spectrum S^(u>) convolved with (T/27r)sinc 2 (wr/27r). 
This is intuitive since one would expect a long r to be 
required to achieve high frequency resolution. In Fig. 
we have applied this to calculate the Mollow spectrum 
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FIG. 4. Comparison of the simulated wave func- 
tion spectrum with the Mollow spectrum convolved with 
(r/27r)sinc 2 (wr/27r) for r = 4/T and Q = 6F The inset 
shows the same comparison for the short time r = 1/T where 
the sidebands cannot be resolved. 

Unraveling the density operator equation to form 
quantum trajectories is of significance for a number of 
reasons. When the system space is large it can be im- 
possible to computationally store and evolve the density 
matrix and one is then forced to use quantum trajectory 
methods which require only the simulation of wave func- 
tions. The more fundamental aspect p[-jlfj| is that the 
system evolution is conditional on the reservoir record 



and when an appropriate measurement basis is used, the 
system may be continuously localized to a region of its 
Hilbert space. The key requirement for this is that the 
dissipative events must provide information about the 
system state. A simple example is in the case of sponta- 
neous emission where imaging the source of the photon 
allows one to localize the position wave function of the 
atom and use a reduced basis set to track the atomic 
motion Jl0| . 

We have presented in this Letter the measurement ba- 
sis which identifies the energy of dissipative events; an 
intrinsically different kind of information to the time do- 
main approaches. There are numerous examples of phys- 
ical systems in which this information is of interest e.g. 
radiative heating in ultracold collisions [[ll) (where the 
frequency of photons is correlated with the internuclear 
separation) and laser cooling jL2j (where as the atomic 
gas cools photons of higher frequency than the driving 
fields are emitted). We have illustrated here the calcula- 
tion of the fluctuation spectrum, although ordered spec- 
tra may also have intrinsic interest for certain problems. 
A final point of note is that the total evolution time may 
be partitioned into intervals of width r and each solved 
using the method we have presented. Then varying r 
one can sweep continuously between use of the time and 
frequency domains. 
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